# This analysis was conducted using R Studio Version 2021.09.1 for Mac and R version 3.5.2 for Mac
rm(list = ls(all = TRUE))
install.packages("pacman")
pacman::p_load(foreign, ggplot2, plm, reshape2, countrycode, sandwich, lmtest, MASS, glmnet, mgcv,
               rworldmap, RColorBrewer, states, mice, VIM, stargazer, margins, clusterSEs, lme4, optimx,
               coefplot, topicmodels, tm, tidytext, dplyr, wordcloud, class, interplot, DataCombine, bife,
               interplot, ivmodel, tidyr, haven, sjPlot, panelView, grid, gridExtra, lattice, zoo, pglm, logistf) # load relevant packages

## Write functions for plotting interactions; used later. Taken from Strezhnev -- https://gist.github.com/fghjorth/45c81a15f40660bc2852
interaction_plot_binary <- function(model, effect, moderator, interaction, varcov="default", conf=.95, title="Marginal effects plot", xlabel="Value of moderator", ylabel="Estimated marginal coefficient", factor_labels=c(0,1)){
  
  # Extract Variance Covariance matrix
  if (varcov == "default"){
    covMat = vcov(model)
  }else{
    covMat = varcov
  }
  
  # Extract the data frame of the model
  mod_frame = model.frame(model)
  
  # Get coefficients of variables
  beta_1 = model$coefficients[[effect]]
  beta_3 = model$coefficients[[interaction]]
  
  # Create list of moderator values at which marginal effect is evaluated
  x_2 <- c(0,1)
  
  # Compute marginal effects
  delta_1 = beta_1 + beta_3*x_2
  
  # Compute variances
  var_1 = covMat[effect,effect] + (x_2^2)*covMat[interaction, interaction] + 2*x_2*covMat[effect, interaction]
  
  # Standard errors
  se_1 = sqrt(var_1)
  
  # Upper and lower confidence bounds
  z_score = qnorm(1 - ((1 - conf)/2))
  upper_bound = delta_1 + z_score*se_1
  lower_bound = delta_1 - z_score*se_1
  
  # Determine the bounds of the graphing area
  max_y = 1
  min_y = 0
  
  # Initialize plotting window
  plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(-.5, 1.5), xlab=xlabel, ylab=ylabel, main=title, xaxt="n")
  
  # Plot points of estimated effects
  points(x=x_2, y=delta_1, pch=16)
  
  # Plot lines of confidence intervals
  lines(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), lty=1)
  points(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), pch=c(25,24), bg="black")
  lines(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), lty=1)
  points(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), pch=c(25,24), bg="black")
  
  # Label the axis
  axis(side=1, at=c(0,1), labels=factor_labels)
  
  # Add a dashed horizontal line for zero
  abline(h=0, lty=3)
  
}
interaction_plot_continuous <- function(model, effect, moderator, interaction, varcov="default", minimum="min", maximum="max", incr="default", num_points = 50, conf=.95, mean=FALSE, median=FALSE, alph=80, rugplot=T, histogram=T, title="Marginal effects plot", xlabel="Value of moderator", ylabel="Estimated marginal coefficient"){
  
  # Define a function to make colors transparent
  makeTransparent<-function(someColor, alpha=alph){
    newColor<-col2rgb(someColor)
    apply(newColor, 2, function(curcoldata){rgb(red=curcoldata[1], green=curcoldata[2],
                                                blue=curcoldata[3],alpha=alpha, maxColorValue=255)})
  }
  
  # Extract Variance Covariance matrix
  if (varcov == "default"){
    covMat = vcov(model)
  }else{
    covMat = varcov
  }
  
  # Extract the data frame of the model
  mod_frame = model.frame(model)
  
  # Get coefficients of variables
  beta_1 = model$coefficients[[effect]]
  beta_3 = model$coefficients[[interaction]]
  
  # Set range of the moderator variable
  # Minimum
  if (minimum == "min"){
    min_val = min(mod_frame[[moderator]])
  }else{
    min_val = minimum
  }
  # Maximum
  if (maximum == "max"){
    max_val = max(mod_frame[[moderator]])
  }else{
    max_val = maximum
  }
  
  # Check if minimum smaller than maximum
  if (min_val > max_val){
    stop("Error: Minimum moderator value greater than maximum value.")
  }
  
  # Determine intervals between values of the moderator
  if (incr == "default"){
    increment = (max_val - min_val)/(num_points - 1)
  }else{
    increment = incr
  }
  
  # Create list of moderator values at which marginal effect is evaluated
  x_2 <- seq(from=min_val, to=max_val, by=increment)
  
  # Compute marginal effects
  delta_1 = beta_1 + beta_3*x_2
  
  # Compute variances
  var_1 = covMat[effect,effect] + (x_2^2)*covMat[interaction, interaction] + 2*x_2*covMat[effect, interaction]
  
  # Standard errors
  se_1 = sqrt(var_1)
  
  # Upper and lower confidence bounds
  z_score = qnorm(1 - ((1 - conf)/2))
  upper_bound = delta_1 + z_score*se_1
  lower_bound = delta_1 - z_score*se_1
  
  # Determine the bounds of the graphing area
  max_y = 1
  min_y = 0
  
  # Make the histogram color
  hist_col = makeTransparent("grey")
  
  # Initialize plotting window
  plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(min_val, max_val), xlab=xlabel, ylab=ylabel, main=title)
  
  # Plot estimated effects
  lines(y=delta_1, x=x_2)
  lines(y=upper_bound, x=x_2, lty=2)
  lines(y=lower_bound, x=x_2, lty=2)
  
  # Add a dashed horizontal line for zero
  abline(h=0, lty=3)
  
  # Add a vertical line at the mean
  if (mean){
    abline(v = mean(mod_frame[[moderator]]), lty=2, col="red")
  }
  
  # Add a vertical line at the median
  if (median){
    abline(v = median(mod_frame[[moderator]]), lty=3, col="blue")
  }
  
  # Add Rug plot
  if (rugplot){
    rug(mod_frame[[moderator]])
  }
  # Add Histogram (Histogram only plots when minimum and maximum are the min/max of the moderator)
  if (histogram & minimum=="min" & maximum=="max"){
    par(new=T)
    hist(mod_frame[[moderator]], axes=F, xlab="", ylab="",main="", border=hist_col, col=hist_col)
  }
}
## Write function to get the mode of data; used later
getmode <- function(v) {
  uniqv <- unique(v)
  uniqv[which.max(tabulate(match(v, uniqv)))]
}
# set working directory (change to your directory)
setwd("ENTER DIRECTORY")

#### Create datasets (ONLY TO REPLICATE DATASET CREATION) ####
## These datasets are loaded later, so this can be skipped for replication of analysis only.
##->## Create IBRD power asymmetry variable
ibrdvote <- read.csv("ibrdvotes.csv") # loads in data on IBRD votes and shares, which was hand-coded from annual reports
# Voting power in the IBRD is revised periodically by the Board of Governor's, and data on voting 
# power can be found in annual reports. Values were hand-coded from these reports
ibrdvote$ctry <- as.character(ibrdvote$ctry) # change format
ibrdvote$ccode <- countrycode(ibrdvote$ctry, "country.name", "iso3c") # add ccodes
gdp <- read.csv("gdp.csv") # load GDP data from WDI
gdp <- gdp[-c(2,3)] # delete extra cols
colnames(gdp) <- c("year", "ccode", "gdp") # rename cols
ibrdgdp <- merge(ibrdvote, gdp, by = c("ccode", "year"), all.x = TRUE) # merge data
ibrdgdp <- ibrdgdp[!ibrdgdp$year %in% 2018:2019,] # remove years for which gdp data systematically not available
ibrdgdpimp <- mice(ibrdgdp, method = "cart", seed = 500, m = 1,
                   pred=quickpred(ibrdgdp, exclude=c("ctry"))) # impute remaining missing gdp data so data is complete
ibrdgdpimp <- mice::complete(ibrdgdpimp)
for (i in 1:7286) {
  year = ibrdgdpimp$year[i]
  temp = sum(ibrdgdpimp$gdp[ibrdgdpimp$year == year])
  ibrdgdpimp$gdpshare[i] <- (ibrdgdpimp$gdp[i] / temp)  
} # calculate gdp share
ibrdgdpimp$asymm <- ibrdgdpimp$gdpshare - ibrdgdpimp$voteshare # calculate asymmetry
# write.csv(ibrdgdpimp, "powerasymmwb.csv") # writes data to file; toggled off for now

##->## Load covariate data and format it appropriately 
wdi <- read.csv("wdi_11.7.19.csv") # load data for IMF performance DV and relevant controls downloaded from WDI on 11/7/19 https://databank.worldbank.org/source/world-development-indicators
wdi <- wdi[-c(2)] # delete extra column
colnames(wdi) <- c("year", "ctry", "ccode", "gdp", "gdpgrow", "reserves", "receipts", "payments", "curraccgdp", "pop") # rename columns
wdi$year <- as.numeric(as.character(wdi$year)) # change data format for variables
wdi$ccode <- as.character(wdi$ccode)
wdi$gdp <- as.numeric(as.character(wdi$gdp))
wdi$gdpgrow <- as.numeric(as.character(wdi$gdpgrow))
wdi$reserves <- as.numeric(as.character(wdi$reserves))
wdi$receipts <- as.numeric(as.character(wdi$receipts))
wdi$payments <- as.numeric(as.character(wdi$payments))
wdi$curraccgdp <- as.numeric(as.character(wdi$curraccgdp))
wdi$pop <- as.numeric(as.character(wdi$pop))
load("Dyadicdata.RData") # load dyadic UN voting data
UNGA <- x # assign to name
UNGA_US <- UNGA[UNGA$ccode1 == 2,] # subset to relations with US
UNGA_US$ccode <- countrycode(UNGA_US$ccode2, origin = "cown", destination = "iso3c") # build ccode
UNGA_US <- UNGA_US[-c(1:2,4:5,7:21)] # delete excess columns
UNGA_US <- aggregate(list(UNGA_US$absidealdiff, UNGA_US$absidealimportantdiff), 
                     by = list(UNGA_US$ccode, UNGA_US$year), FUN = mean) # eliminates issue with duplicates
colnames(UNGA_US) <- c("ccode", "year", "absidealdiff", "absidealimportantdiff") # rename columns
wdi <- merge(wdi, UNGA_US, by = c("ccode","year"), all.x = T) # merge
polity <- read.csv("p4v2017.csv") # load polity data
polity$ccode <- countrycode(polity$ctry, origin = "country.name", destination = "iso3c") # build ccode
polity <- polity[-c(2)] # delete unnecessary column
polity <- aggregate(list(polity$polity2), 
                    by = list(polity$ccode, polity$year), FUN = mean) # eliminates issue with duplicates
colnames(polity) <- c("ccode", "year", "polity2") # rename columns
wdi <- merge(wdi, unique(polity), by = c("ccode","year"), all.x = T) # merge  
usaid <- read.csv("us_foreignaid_greenbook.csv") # load foreign aid greenbook data downloaded August 2020
usaid <- usaid[-c(2,4:8)] # delete extra cols
colnames(usaid) <- c("year", "ctry", "USaid") # rename cols
usaid <- aggregate(usaid$USaid, by = list(usaid$year, usaid$ctry), FUN = sum) # aggregate by country
colnames(usaid) <- c("year", "ctry", "USaid") # change col names
usaid$ccode <- countrycode(usaid$ctry, origin = "country.name", destination = "iso3c") # add country code
usaid <- usaid[-c(2)] # delete col
wdi <- merge(wdi, unique(usaid), by = c("ccode", "year"), all.x = TRUE) # merge data
ustrade <- read.csv("Dyadic_COW_4.0.csv") # load COW 4.0 dyadic trade data downloaded January 2020 https://correlatesofwar.org/data-sets/bilateral-trade
ustrade <- ustrade[ustrade$ccode1 == 2,] # restrict to dyadic with US
ustrade$tottrade <- ustrade$flow1 + ustrade$flow2 # calculate total trade
ustrade <- ustrade[c(2,3,24)] # delete extra cols
colnames(ustrade) <- c("ccode", "year", "tottrade") # rename columns
ustrade$ccode <- countrycode(ustrade$ccode, origin = "cown", destination = "iso3c") # generae ccode
wdi <- merge(wdi, unique(ustrade), by = c("ccode", "year"), all.x = TRUE) # merge data
allyfull <- read.csv("ally.csv") # load ATOP alliance data downloaded in December 2019 http://www.atopdata.org/
allyfull <- allyfull[allyfull$ccode1 == 2,] # subset to only allies of the US
allyfull <- allyfull[-c(1:3, 5:13,19)] # subset to relevant columns
colnames(allyfull)[1] <- "ccode" # change col name
allyfull <- aggregate(list(allyfull$defense, allyfull$neutrality, allyfull$nonaggression, allyfull$entente), by = 
                        list(allyfull$ccode, allyfull$year), FUN = max) # eliminate duplicates
colnames(allyfull) <- c("ccode", "year", "defense", "neutrality", "nonaggression", "entente") # change col names
allyfull$ccode <- countrycode(allyfull$ccode, origin = "cown", destination = "iso3c") # harmonize ccodes
wdi <- merge(wdi, allyfull, by = c("ccode", "year"), all.x = T) # merge
wdi$defense[is.na(wdi$defense)] <- 0 # NAs represent a lack of alliance, so recode for all alliance types
wdi$neutrality[is.na(wdi$neutrality)] <- 0
wdi$entente[is.na(wdi$entente)] <- 0
wdi$nonaggression[is.na(wdi$nonaggression)] <- 0
unsc <- read.csv("unsc.csv") # load data on UNSC membership from Dreher, Sturm, and Vreeland (2009) https://www.uni-heidelberg.de/fakultaeten/wiso/awi/professuren/intwipol/datasets_en.html
# downloaded in November 2018; I hand-coded the more recent years from UN online materials to bring the data up to date
unsc$ccode <- countrycode(unsc$ctry, origin = "country.name", destination = "iso3c") # harmonize country codes
unsc$unsc[unsc$unsc == "."] <- 0 # change empty observations to 0
unsc <- na.omit(unsc) # omit NAs (states that no longer exist)
unsc <- unsc[-c(1)] # delete extra col
unsc$unsc <- as.numeric(as.character(unsc$unsc)) # change data format
unsc <- aggregate(unsc$unsc, by = list(unsc$ccode, unsc$year), FUN = max) # eliminate issues with duplicates
colnames(unsc) <- c("ccode", "year", "unsc") # change col names
wdi <- merge(unique(wdi), unique(unsc), by = c("ccode", "year"), all.x = TRUE) # merge data
wdi$unsc[is.na(wdi$unsc)] <- 0 # change NAs to 0 because missingness denotes not a member
wdi <- wdi[-c(1:5),] # deletes blank rows
ids <- read.csv("IDS.csv") # load international debt statistics downloaded from World Bank in June 2020
ids <- ids[-c(2,3)] # delete extra column
colnames(ids) <- c("year", "ccode", "ibrdcom", "privcom", "offcom")
ids$totcom <- as.numeric(ids$privcom) + as.numeric(ids$offcom)
ids$year <- as.numeric(ids$year) # change data format
ids$privcom <- as.numeric(ids$privcom) # change data format
ids$ibrdcom <- as.numeric(ids$ibrdcom) # change data format
ids$offcom <- as.numeric(ids$offcom) # change data format
ids[is.na(ids)] <- 0 # if not in data, commitments are 0
wdi <- merge(unique(wdi), unique(ids), by = c("ccode", "year"), all.x = TRUE)
# write.csv(wdi, "covariates.csv") # writes to file; toggled off for now
wdiimp <- mice(wdi, method = "cart", seed = 500, m = 1,
                  pred=quickpred(wdi, exclude=c("ctry", "ccode"))) # impute missing data
wdiimp <- mice::complete(wdiimp) 
# write.csv(wdiimp, "covariatesimp.csv") # writes to file; toggled off for now

##->## Load additional datasets needed for analysis
refp <- read.csv("reformp.csv") # load hand-coded IO-year data on policy reforms
ibrdvotes <- read.csv("ibrdvotes.csv") # load hand-coded ctry-year data on ibrd country-specific reforms
wbasymm <- read.csv("powerasymmwb.csv") # load ibrd power asymmetry data, which was created above
devmem <- read.csv("development_membership.csv") # load hand-coded data on membership in development IOs (many are missing from Pevehouse data so I prefer this method)
perfwbrate <- read.csv("IEG_clean.csv") # load WB performance rating data

##->## Merge together various datasets for IBRD analysis
devmem <- devmem[-c(2,3)] # delete extra cols
wbasymm <- wbasymm[-c(1,4,7,8)] # delete extra cols
devmem$ccode <- countrycode(devmem$ccode, origin = "cown", destination = "iso3c") # harmonize ccodes
devmemvote <- merge(unique(wbasymm), unique(devmem), by = c("ccode", "year")) # merge data
devmemvote <- devmemvote[!is.na(devmemvote$ccode),] # eliminate NAs
devmemvote <- unique(devmemvote) # eliminate any duplicates
devmemvotefull <- merge(unique(devmemvote), unique(wdiimp), by = c("ccode", "year"))
tempdv <- devmemvotefull[c(1:4)] # create temp dataframe for lagging
devmemvotefulllag <- devmemvotefull # lag data by one year
devmemvotefulllag$year <- devmemvotefulllag$year + 1
devmemvotefulllag <- devmemvotefulllag[-c(3,4)] # delete extra cols
devmemvotefulllag <- merge(devmemvotefulllag, tempdv, by = c("ccode", "year")) # merge data
ibrdvotes$ccode <- countrycode(ibrdvotes$ctry, origin = "country.name", destination = "cown") # harmonize ccodes
ibrdvotes$polrefben[is.na(ibrdvotes$polrefben)] <- 0 # change NAs to 0 because they are true zeroes
ctryref <- ibrdvotes[-c(1,3:4)] # delete extra cols
ctryref$ccode <- countrycode(ctryref$ccode, origin = "cown", destination = "iso3c")
devmemvotefulllag3 <- merge(devmemvotefulllag, ctryref, by = c("ccode", "year")) # merge data
devmemvotefulllag3$ccode <- as.character(devmemvotefulllag3$ccode) # change format
devmemvotefulllag3$ccode <- countrycode(devmemvotefulllag3$ccode, origin = "iso3c", destination = "cown") # change format
devmemvotefulllag3$ccode <- as.numeric(devmemvotefulllag3$ccode)
attach(devmemvotefulllag3) # attach frame to create reform beneficiary data
devmemvotefulllag3.1 <- devmemvotefulllag3 %>% arrange(ccode, year) # order data
devmemvotefulllag3.1 <- change(data = devmemvotefulllag3.1, Var = "voteshare", GroupVar = "ccode", NewVar = "votech", 
                               slideBy = -1, type = "absolute") # code countries with positive vote share changes
devmemvotefulllag3.1$votechpos <- 0
devmemvotefulllag3.1$votechpos[devmemvotefulllag3.1$votech > 0] <- 1
devmemvotefulllag3.1$polrefbenall <- 0 # code countries with positive reforms OR vote share changes
devmemvotefulllag3.1$polrefbenall <- devmemvotefulllag3.1$votechpos + devmemvotefulllag3.1$polrefben
devmemvotefulllag3.1$polrefbenall[devmemvotefulllag3.1$polrefbenall > 1] <- 1
devmemvotefulllag3.1$polrefben[devmemvotefulllag3.1$year %in% c(1959, 1970, 1984, 1986, 2010, 2011, 
                                                                2018) & devmemvotefulllag3.1$votech > 0] <- 1 # code countries with positive reforms OR vote share changes only for ambiguous reform years
perfwbrate <- perfwbrate[perfwbrate$io=="IBRD",] # subset performance data to IBRD projects
perfwbrate$ccode <- countrycode(perfwbrate$ctry, origin = "country.name", destination = "cown") # harmonize ccodes
perfwbrateagg <- aggregate(perfwbrate$rate_num, by = list(perfwbrate$ccode, perfwbrate$year), FUN = mean) # aggregate performance by ctry-year
colnames(perfwbrateagg) <- c("ccode", "year", "ratenum") # change colnames
perfwbrateaggyr <- aggregate(perfwbrate$rate_num, by = list(perfwbrate$year), FUN = mean) # aggregate performance by year
colnames(perfwbrateaggyr) <- c("year", "ratenum") # change colnames
perfwbrateaggyr$year <- as.numeric(perfwbrateaggyr$year) # change format
perfwbrate1 <- perfwbrateagg # duplicate for lagged data
perfwbrate1$year <- perfwbrate1$year + 1 # add year for lagged data
ratefull1 <- merge(devmemvotefulllag3.1, perfwbrate1, by = c("ccode", "year"), all =T) # merge data
ratefull1  <- unique(ratefull1) # eliminate any duplicates
# write.csv(ratefull1, "ibrdreformfull.csv") # write to file; turned off for now

##->## Do the same without imputed covariates
devmemvotefull <- merge(devmemvote, wdi, by = c("ccode", "year"))
tempdv <- devmemvotefull[c(1:4)] # create temp dataframe for lagging
devmemvotefulllag <- devmemvotefull # lag data by one year
devmemvotefulllag$year <- devmemvotefulllag$year + 1
devmemvotefulllag <- devmemvotefulllag[-c(3,4)] # delete extra cols
devmemvotefulllag <- merge(devmemvotefulllag, tempdv, by = c("ccode", "year")) # merge data
ibrdvotes$ccode <- countrycode(ibrdvotes$ctry, origin = "country.name", destination = "cown") # harmonize ccodes
ibrdvotes$polrefben[is.na(ibrdvotes$polrefben)] <- 0 # change NAs to 0 because they are true zeroes
ctryref <- ibrdvotes[-c(1,3:4)] # delete extra cols
ctryref$ccode <- countrycode(ctryref$ccode, origin = "cown", destination = "iso3c") # harmonize ccodes
devmemvotefulllag3 <- merge(devmemvotefulllag, ctryref, by = c("ccode", "year")) # merge data
devmemvotefulllag3$ccode <- as.character(devmemvotefulllag3$ccode) # change format
devmemvotefulllag3$ccode <- countrycode(devmemvotefulllag3$ccode, origin = "iso3c", destination = "cown") # change format
devmemvotefulllag3$ccode <- as.numeric(devmemvotefulllag3$ccode) # change format
attach(devmemvotefulllag3) # attach frame to create reform beneficiary data
devmemvotefulllag3.1 <- devmemvotefulllag3 %>% arrange(ccode, year) # order data
devmemvotefulllag3.1 <- change(data = devmemvotefulllag3.1, Var = "voteshare", GroupVar = "ccode", NewVar = "votech", 
                               slideBy = -1, type = "absolute") # code countries with positive vote share changes
devmemvotefulllag3.1$votechpos <- 0 # create variable for vote change position
devmemvotefulllag3.1$votechpos[devmemvotefulllag3.1$votech > 0] <- 1 # code as 1 for positive change
devmemvotefulllag3.1$polrefbenall <- 0 # code countries with positive reforms OR vote share changes
devmemvotefulllag3.1$polrefbenall <- devmemvotefulllag3.1$votechpos + devmemvotefulllag3.1$polrefben # sum vote change position and beneficial reform variables
devmemvotefulllag3.1$polrefbenall[devmemvotefulllag3.1$polrefbenall > 1] <- 1 # recode as binary
devmemvotefulllag3.1$polrefben[devmemvotefulllag3.1$year %in% c(1959, 1970, 1984, 1986, 2010, 2011, 
                                                                2018) & devmemvotefulllag3.1$votech > 0] <- 1 # code countries with positive reforms OR vote share changes only for ambiguous reform years
ratefull2 <- merge(devmemvotefulllag3.1, perfwbrate1, by = c("ccode", "year"), all =T) # merge data
ratefull2 <- unique(ratefull2) # eliminate any duplicates
# write.csv(ratefull2, "ibrdreformsub.csv") # write to file; turned off for now

#### Descriptive plots ####
##->## Reform frequency over time
refpwb <- refp[refp$io == "IBRD",]
refpwb %>%
  ggplot(aes(x = year, y = reformp)) +
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_line(stat = "identity", position = "dodge") +
  labs(x = "Year", y = "Number of Reforms") +
  scale_x_continuous(breaks = c(1950, 1960, 1970, 1980, 1990, 2000, 2010, 2020)) +
  scale_y_continuous(breaks = c(2, 4, 6, 8, 10 , 12))# Exported and presented as Figure A2

##->## IBRD map of total beneficial reforms 
map2 <- as.data.frame(cbind(as.character(ibrdimp$ccode), as.numeric(as.character(ibrdimp$year)), as.numeric(as.character(ibrdimp$polrefben))))
colnames(map2) <- c("ccode", "year", "reforms")
map2$ctry <- countrycode(map2$ccode, origin = "cown", destination = "country.name")
map2$reforms[is.na(map2$reforms)] <- 0
map2$reforms <- as.numeric(as.character(map2$reforms))
map2$ccode <- as.numeric(as.character(map2$ccode))
map2$year <- as.numeric(as.character(map2$year))
map2 <- aggregate(list(map2$reforms), by = list(map2$ccode, map2$ctry), sum)
colnames(map2) <- c("ccode", "ctry", "reforms")
map2$ccode <- countrycode(map2$ccode, origin = "cown", destination = "iso3c")
maps2 <- joinCountryData2Map(map2, joinCode = "ISO3", nameJoinColumn = "ccode", nameCountryColumn = "ctry")
maps2 <- subset(maps2, continent != "Antarctica")
mapss2 <- mapCountryData(mapToPlot = maps2, nameColumnToPlot = "reforms", mapRegion = "world", 
                         catMethod = c(0:4), colourPalette = brewer.pal(11, name = "Blues"),
                         addLegend = FALSE, mapTitle = "")
do.call( addMapLegend, c(mapss2, legendWidth=0.5, legendMar = 5)) # Exported and presented as Figure 1

##->## Categories of reforms at IMF and WB (original descriptive analysis)
cat <- read.csv("reformp_cat_ibrd_alt.csv") # load hand-coded data on categories of reforms
catagg <- aggregate(list(cat$miss, cat$rule, cat$staff, cat$lend, 
                         cat$oth, cat$vote2, cat$subscr, cat$capstock), by = list(cat$io), FUN = sum)
colnames(catagg) <- c("IO", "Expansion", "Rules/By-Laws", "New Staff", "Lending Policy", 
                      "Other", "Vote Share", "Subscriptions", "Capital Stock") # used this data to hand code categories in more user friendly way, loaded below
catagg2 <- read.csv("catagg2.csv") # load more user friendly data
positions <- c("Lending Policy", "Expansion", "New Staff", "Rules/By-Laws", 
               "Vote Share", "Subscriptions", "Capital Stock", "Other")

ggplot(catagg2, aes(y=freq, x=cat)) +
  theme_bw() + 
  geom_bar(position="dodge", stat="identity") +
  labs(x = "Reform Category", y = "Frequency") +
  scale_x_discrete(limits = positions) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) # Exported and presented as Figure A1

catagg2 <- aggregate(list(cat$miss, cat$rule, cat$staff, cat$lend, 
                         cat$oth, cat$vote), by = list(cat$io, cat$year), FUN = sum)
colnames(catagg2) <- c("IO", "Year", "Expansion", "Rules/By-Laws", "New Staff", "Lending Policy", 
                      "Other", "Vote Share") # used this data to hand code categories in more user friendly way, loaded below

p1 <- ggplot(catagg2, aes(y=Expansion, x=Year)) +
  theme_bw() + 
  geom_line() +
  labs(x = "Year", y = "Frequency", title = "Expansion") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  scale_x_continuous(breaks = c(1950,1960,1970,1980,1990,2000,2010,2020)) +
  scale_y_continuous(breaks=c(1,2))
p2 <- ggplot(catagg2, aes(y=`Rules/By-Laws`, x=Year)) +
  theme_bw() + 
  geom_line() +
  labs(x = "Year", y = "Frequency", title="Rules/By-Laws") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  scale_x_continuous(breaks = c(1950,1960,1970,1980,1990,2000,2010,2020)) +
  scale_y_continuous(breaks=c(1,2))
p3 <- ggplot(catagg2, aes(y=`New Staff`, x=Year)) +
  theme_bw() + 
  geom_line() +
  labs(x = "Year", y = "Frequency", title="New Staff") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  scale_x_continuous(breaks = c(1950,1960,1970,1980,1990,2000,2010,2020)) +
  scale_y_continuous(breaks=c(1,2))
p4 <- ggplot(catagg2, aes(y=`Lending Policy`, x=Year)) +
  theme_bw() + 
  geom_line() +
  labs(x = "Year", y = "Frequency", title="Lending Policy") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  scale_x_continuous(breaks = c(1950,1960,1970,1980,1990,2000,2010,2020)) +
  scale_y_continuous(breaks=c(1,2))
p5 <- ggplot(catagg2, aes(y=`Vote Share`, x=Year)) +
  theme_bw() + 
  geom_line() +
  labs(x = "Year", y = "Frequency", title="Vote Share") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  scale_x_continuous(breaks = c(1950,1960,1970,1980,1990,2000,2010,2020)) +
  scale_y_continuous(breaks=c(2,4,6,8,10))
p6 <- ggplot(catagg2, aes(y=`Other`, x=Year)) +
  theme_bw() + 
  geom_line() +
  labs(x = "Year", y = "Frequency", title="Other") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  scale_x_continuous(breaks = c(1950,1960,1970,1980,1990,2000,2010,2020)) +
  scale_y_continuous(breaks = c(1,2))

grid.arrange(
  p1,
  p2,
  p3,
  p4,
  p5,
  p6,
  nrow = 3
) # Exported and presented as Figure A3

##->## China and India performance over time (for case)

allymean <- ibrdimp[c(3,48,60)]
allymean <- na.omit(allymean)
allymean <- allymean %>%
  group_by(year) %>%
  summarise_at(vars(-defense), funs(mean(.)))
allymean$year <- allymean$year - 1
allymean <- allymean[allymean$year > 1988,] # restrict to period where China also receives loans

chin <- ibrdimp[ibrdimp$ccode == 710 | ibrdimp$ccode == 750,]
chin$ccode[chin$ccode == 710] <- "China"
chin$ccode[chin$ccode == 750] <- "India"
chin$year <- chin$year - 1
chin <- chin[chin$year > 1988,] # restrict to period where China also receives loans
ggplot(chin, aes(y = ratenum, x = year, col = factor(ccode))) +
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  geom_smooth(method = "loess", se = F) +
  labs(x = "Year", y = "Performance", col = "Country") + 
  # geom_line() +
  geom_smooth(data = allymean, aes(y = ratenum, x = year, col = "Global Mean"), se = F) # plots when performance had dipped for these countries relative to global mean; exported as Figure 3

chin2 <- ibrdimp[c(2,3,56)]
chin2 <- chin2[chin2$ccode == 710 | chin2$ccode == 750,] # tells us when China and India received reforms

braz <- ibrdimp[ibrdimp$ccode == 140,] # repeat for brazil; Figure A5
braz$ccode <- "Brazil"
braz$year <- braz$year - 1
braz <- braz[braz$year < 2017,]

allymean <- ibrdimp[c(3,48,63)]
allymean <- na.omit(allymean)
allymean <- allymean %>%
  group_by(year) %>%
  summarise_at(vars(-defense), funs(mean(.)))
allymean$year <- allymean$year - 1

ggplot(braz, aes(y = ratenum, x = year, col = factor(ccode))) +
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.title.x = element_text(size = 22), axis.title.y = element_text(size = 22),
        axis.text.x = element_text(size = 20), axis.text.y = element_text(size = 20),
        legend.text = element_text(size = 20), legend.title = element_text(size = 22)) + 
  geom_smooth(method = "loess", se = F) +
  labs(x = "Year", y = "Performance", col = "Country") + 
  # geom_line() +
  geom_smooth(data = allymean, aes(y = ratenum, x = year, col = "Global Mean"), se = F) # 

# Repeat for Indonesia; Figure A6

indo <- ibrdimp[ibrdimp$ccode == 850,]
indo$ccode <- "Indonesia"
indo$year <- indo$year - 1
indo <- indo[indo$year < 2017,]

allymean <- ibrdimp[c(3,48,63)]
allymean <- na.omit(allymean)
allymean <- allymean %>%
  group_by(year) %>%
  summarise_at(vars(-defense), funs(mean(.)))
allymean$year <- allymean$year - 1

ggplot(indo, aes(y = ratenum, x = year, col = factor(ccode))) +
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.title.x = element_text(size = 22), axis.title.y = element_text(size = 22),
        axis.text.x = element_text(size = 20), axis.text.y = element_text(size = 20),
        legend.text = element_text(size = 20), legend.title = element_text(size = 22)) + 
  geom_smooth(method = "loess", se = F) +
  labs(x = "Year", y = "Performance", col = "Country") + 
  # geom_line() +
  geom_smooth(data = allymean, aes(y = ratenum, x = year, col = "Global Mean"), se = F) # 

# Repeat for Eastern bloc; Figure A7

east <- ibrdimp[ibrdimp$ccode == 371 | ibrdimp$ccode == 373 | ibrdimp$ccode == 370 | ibrdimp$ccode == 366 | 
                  ibrdimp$ccode == 372 | ibrdimp$ccode == 705 | ibrdimp$ccode == 703 | ibrdimp$ccode == 367 |
                  ibrdimp$ccode == 368 | ibrdimp$ccode == 359 | ibrdimp$ccode == 365 | ibrdimp$ccode == 702 |
                  ibrdimp$ccode == 701 | ibrdimp$ccode == 369 | ibrdimp$ccode == 704,]
east$ccode[east$ccode == 371] <- "Armenia"
east$ccode[east$ccode == 373] <- "Azerbaijan"
east$ccode[east$ccode == 370] <- "Belarus"
east$ccode[east$ccode == 366] <- "Estonia"
east$ccode[east$ccode == 372] <- "Georgia"
east$ccode[east$ccode == 705] <- "Kazakhstan"
east$ccode[east$ccode == 703] <- "Kyrgyzstan"
east$ccode[east$ccode == 367] <- "Latvia"
east$ccode[east$ccode == 368] <- "Lithuania"
east$ccode[east$ccode == 359] <- "Moldova"
east$ccode[east$ccode == 365] <- "Russia"
east$ccode[east$ccode == 702] <- "Tajikistan"
east$ccode[east$ccode == 701] <- "Turkmenistan"
east$ccode[east$ccode == 369] <- "Ukraine"
east$ccode[east$ccode == 704] <- "Uzbekistan"
east$year <- east$year - 1
eastagg <- aggregate(east$ratenum, by = list(east$year), FUN = mean, na.rm=T)
colnames(eastagg) <- c("year", "ratenum")
ggplot(eastagg, aes(y = ratenum, x = year, col = "Bloc Mean")) +
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.title.x = element_text(size = 22), axis.title.y = element_text(size = 22),
        axis.text.x = element_text(size = 20), axis.text.y = element_text(size = 20),
        legend.text = element_text(size = 20), legend.title = element_text(size = 22)) + 
  geom_smooth(method = "loess", se = F) +
  labs(x = "Year", y = "Performance", col = "Country") + 
  scale_x_continuous(limits = c(1995, 2020)) +
  geom_smooth(data = allymean, aes(y = ratenum, x = year, col = "Global Mean"), se = F) # plots when performance had dipped for these countries relative to global mean; exported as Figure 5

# Raw associations between ally, perf, and reform; see p. 23 for discussion
cor.test(ibrdimp$defense, ibrdimp$polrefben, method = "pearson")
mean(ibrdimp$polrefben[ibrdimp$defense == 1])
mean(ibrdimp$polrefben[ibrdimp$defense == 0])
quantile(ibrdimp$ratenum[ibrdimp$defense == 0], na.rm=T)
q <- ibrdimp %>% 
  group_by(year) %>% 
  summarise(value.25th = quantile(ratenum, 0.25, na.rm=T))
ibrdtemp <- merge(ibrdimp, q, by = "year")
ibrdtemp <- ibrdtemp[ibrdtemp$nodefense == 1,]
mean(ibrdtemp$polrefben[ibrdtemp$ratenum <= ibrdtemp$value.25th], na.rm=T)
mean(ibrdtemp$polrefben[ibrdtemp$ratenum > ibrdtemp$value.25th], na.rm=T)

#### Data pre-processing (START HERE FOR ANALYSIS REPLICATION) ####
ibrdimp <- read.csv("ibrdreformfull.csv") # load in the datasets which were built above
ibrdnoimp <- read.csv("ibrdreformsub.csv")

ibrdimp <- ibrdimp[!is.na(ibrdimp$polrefben),]
ibrdimp$gdppc <- log1p(ibrdimp$gdp / ibrdimp$pop)
ibrdimp$gdp <- log1p(ibrdimp$gdp / 1000000) # set negatives to, divide values by one million, and log
ibrdimp$reserves[ibrdimp$reserves < 0] <- 0
ibrdimp$reserves <- log1p(ibrdimp$reserves / 1000000)
ibrdimp$payments[ibrdimp$payments < 0] <- 0
ibrdimp$payments <- log1p(ibrdimp$payments / 1000000)
ibrdimp$receipts[ibrdimp$receipts < 0] <- 0
ibrdimp$receipts <- log1p(ibrdimp$receipts / 1000000)
ibrdimp$ibrdcom <- log1p(ibrdimp$ibrdcom / 1000000)
ibrdimp$totcom <- log1p(ibrdimp$totcom / 1000000)
ibrdimp$pop <- log1p(ibrdimp$pop / 1000000)
ibrdimp$USaid[ibrdimp$USaid < 0] <- 0
ibrdimp$USaid <- log1p(ibrdimp$USaid)
ibrdimp$USaid[is.na(ibrdimp$USaid)] <- 0
ibrdimp$tottrade[ibrdimp$tottrade < 0] <- 0
ibrdimp$tottrade <- log1p(ibrdimp$tottrade)
ibrdimp$nodefense <- ifelse(ibrdimp$defense == 1, 0, 1)
colnames(ibrdimp)[34] <- "numIOs"
defagg <- ibrdimp %>% 
  group_by(ccode) %>% 
  summarise(defmode = getmode(defense)) # calculate whether ally or not is the modal position by country
ibrdimp <- left_join(ibrdimp, defagg, by = "ccode")
ibrdimp$nodefmode <- ifelse(ibrdimp$defmode == 1, 0, 1)
neutagg <- ibrdimp %>% 
  group_by(ccode) %>% 
  summarise(neutmode = getmode(neutrality)) # calculate whether ally or not is the modal position by country
ibrdimp <- left_join(ibrdimp, neutagg, by = "ccode")
ibrdimp$noneutmode <- ifelse(ibrdimp$neutmode == 1, 0, 1)
ententeagg <- ibrdimp %>% 
  group_by(ccode) %>% 
  summarise(ententemode = getmode(entente)) # calculate whether ally or not is the modal position by country
ibrdimp <- left_join(ibrdimp, ententeagg, by = "ccode")
ibrdimp$noententemode <- ifelse(ibrdimp$ententemode == 1, 0, 1)
ibrdimp$noentente <- ifelse(ibrdimp$entente == 1, 0, 1)
nonaggagg <- ibrdimp %>% 
  group_by(ccode) %>% 
  summarise(nonaggmode = getmode(nonaggression)) # calculate whether ally or not is the modal position by country
ibrdimp <- left_join(ibrdimp, nonaggagg, by = "ccode")
ibrdimp$nononaggmode <- ifelse(ibrdimp$nonaggmode == 1, 0, 1)
ibrdimp$nononagg <- ifelse(ibrdimp$nonaggression == 1, 0, 1)
ibrdimp <- ibrdimp[ibrdimp$ccode != 2,]

ibrdnoimp <- ibrdnoimp[!is.na(ibrdnoimp$polrefben),] # repeat for other data sets
ibrdnoimp$gdppc <- log1p(ibrdnoimp$gdp / ibrdnoimp$pop)
ibrdnoimp$gdp <- log1p(ibrdnoimp$gdp / 1000000)
ibrdnoimp$reserves[ibrdnoimp$reserves < 0] <- 0
ibrdnoimp$reserves <- log1p(ibrdnoimp$reserves / 1000000)
ibrdnoimp$payments[ibrdnoimp$payments < 0] <- 0
ibrdnoimp$payments <- log1p(ibrdnoimp$payments / 1000000)
ibrdnoimp$receipts[ibrdnoimp$receipts < 0] <- 0
ibrdnoimp$receipts <- log1p(ibrdnoimp$receipts / 1000000)
ibrdnoimp$ibrdcom <- log1p(ibrdnoimp$ibrdcom / 1000000)
ibrdnoimp$totcom <- log1p(ibrdnoimp$totcom / 1000000)
ibrdnoimp$pop <- log1p(ibrdnoimp$pop / 1000000)
ibrdnoimp$USaid[ibrdnoimp$USaid < 0] <- 0
ibrdnoimp$USaid <- log1p(ibrdnoimp$USaid)
ibrdnoimp$USaid[is.na(ibrdnoimp$USaid)] <- 0
ibrdnoimp$tottrade[ibrdnoimp$tottrade < 0] <- 0
ibrdnoimp$tottrade <- log1p(ibrdnoimp$tottrade)
ibrdnoimp$nodefense <- ifelse(ibrdnoimp$defense == 1, 0, 1)
colnames(ibrdnoimp)[34] <- "numIOs"
defagg <- ibrdnoimp %>% 
  group_by(ccode) %>% 
  summarise(defmode = getmode(defense)) # calculate whether ally or not is the modal position by country
ibrdnoimp <- left_join(ibrdnoimp, defagg, by = "ccode")
ibrdnoimp$nodefmode <- ifelse(ibrdnoimp$defmode == 1, 0, 1)
ibrdnoimp <- ibrdnoimp[ibrdnoimp$ccode != 2,]
ibrdnoimp$offcom[is.na(ibrdnoimp$offcom)] <- 0
ibrdnoimp$totcom[is.na(ibrdnoimp$totcom)] <- 0
ibrdnoimp$privcom[is.na(ibrdnoimp$privcom)] <- 0
ibrdnoimp$ibrdcom[is.na(ibrdnoimp$ibrdcom)] <- 0

#### Descriptive stats ####
desc <- read.csv("ibrdreformfull.csv") # load in the datasets which were built above
desc <- desc[!is.na(desc$polrefben),]
desc <- desc[-c(1:4,6:33,35,36,41,43,49:51,54,55,57,58,60:62)]
desc$reserves <- desc$reserves / 1000000
desc$reserves[desc$reserves < 0] <- 0
desc$receipts <- desc$receipts / 1000000
desc$receipts[desc$receipts < 0] <- 0
desc$payments <- desc$payments / 1000000
desc$payments[desc$payments < 0] <- 0
desc$pop <- desc$pop / 1000000
desc$USaid <- desc$USaid / 1000000
desc$USaid[desc$USaid < 0] <- 0
desc$ibrdcom <- desc$ibrdcom / 1000000
desc$totcom <- desc$totcom / 1000000
desc$tottrade[desc$tottrade < 0] <- 0
covs_desc <- c("Vote-power asymmetry", "Number IOs", "GDP growth", "Reserves", "Receipts", "Payments", "Population",
               "UN voting (ideal pt dist from U.S.)", "Polity2", "U.S. aid", "Total trade", "Ally", "UNSC member",
               "World Bank commitments", "Bilateral and private commitments", "Number of reforms", "Performance")
stargazer(desc, omit.summary.stat = c("p25", "p75"), digits = 2, style = "ajps",
          covariate.labels = covs_desc) # print table; presented as Table A3

desc <- read.csv("ibrdreformsub.csv") # load in the datasets which were built above
desc <- desc[!is.na(desc$polrefben),]
desc <- desc[-c(1:4,6:33,35,36,41,43,49:51,54,55,57,58,60:62)]
desc$reserves <- desc$reserves / 1000000
desc$reserves[desc$reserves < 0] <- 0
desc$receipts <- desc$receipts / 1000000
desc$receipts[desc$receipts < 0] <- 0
desc$payments <- desc$payments / 1000000
desc$payments[desc$payments < 0] <- 0
desc$pop <- desc$pop / 1000000
desc$USaid <- desc$USaid / 1000000
desc$USaid[desc$USaid < 0] <- 0
desc$ibrdcom <- desc$ibrdcom / 1000000
desc$ibrdcom[is.na(desc$ibrdcom)] <- 0
desc$totcom <- desc$totcom / 1000000
desc$totcom[is.na(desc$totcom)] <- 0
desc$tottrade[desc$tottrade < 0] <- 0
covs_desc <- c("Vote-power asymmetry", "Number IOs", "GDP growth", "Reserves", "Receipts", "Payments", "Population",
               "UN voting (ideal pt dist from U.S.)", "Polity2", "U.S. aid", "Total trade", "Ally", "UNSC member",
               "World Bank commitments", "Bilateral and private commitments", "Number of reforms", "Performance")
stargazer(desc, omit.summary.stat = c("p25", "p75"), digits = 2, style = "ajps",
          covariate.labels = covs_desc) # print table; presented as Table A2

pr <- mean(abs(ibrdimp$votech[(ibrdimp$year == 2011 | ibrdimp$year == 2012) & ibrdimp$votech > 0]), na.rm=T)
po <- mean(abs(ibrdimp$votech[ibrdimp$year < 2011]), na.rm=T)

pr/po # 2010 reforms were 8% larger than those preceding it; discussed in paper as watershed

#### Main Analyses ####
# Bivariate tests with 1 year lag
ratebivlag1 <- plm(polrefben ~ nodefense*ratenum, data = ibrdnoimp, 
                   index = c("ccode", "year"), model = "within", effect = "twoways")
(seratebivlag1 <- coeftest(ratebivlag1, vcov=vcovHC(ratebivlag1, cluster="group"))) # interaction model

ratebivlag2 <- plm(polrefben ~ defense, data = ibrdnoimp, 
                   index = c("ccode", "year"), model = "within", effect = "twoways")
(seratebivlag2 <- coeftest(ratebivlag2, vcov=vcovHC(ratebivlag2, cluster="group"))) # ally only model

rpss <- plm(polrefben ~ nodefense + ratenum, data = ibrdimp, 
            index = c("ccode", "year"), model = "within",
            effect = "twoways") # simple model without interaction for t test of interaction term
(serpss <- coeftest(rpss, vcov=vcovHC(rpss, cluster="group"))) # ally only model

waldtest(ratebivlag1, rpss, vcov = vcovHC) # p value is 0.016

# Partial cov tests with 1 year lag
ratepartlag1 <- plm(polrefben ~ nodefense*ratenum + numIOs + asymm + unsc +
                      polity2 + absidealimportantdiff + USaid + tottrade, data = ibrdimp, 
                    index = c("ccode", "year"), model = "within",
                    effect = "twoways")
(seratepartlag1 <- coeftest(ratepartlag1, vcov=vcovHC(ratepartlag1, cluster="group"))) # interaction model

interaction_plot_continuous(model = ratepartlag1, conf = 0.9, effect = "nodefense", moderator = "ratenum", 
                            interaction = "nodefense:ratenum", title = "", xlabel = "Performance", 
                            ylabel = "Probability of reform for non-ally",
                            varcov = vcovHC(ratepartlag1, cluster="group")) #  exported and presented as Figure 4

ratepartlag2 <- plm(polrefben ~ defense + numIOs + asymm + unsc +
                      polity2 + absidealimportantdiff + USaid + tottrade, data = ibrdimp, 
                    index = c("ccode", "year"), model = "within",
                    effect = "twoways")
(seratepartlag2 <- coeftest(ratepartlag2, vcov=vcovHC(ratepartlag2, cluster="group"))) # ally only model

# Full tests with 1 year lag
ratefulllag1 <- plm(polrefben ~ nodefense*ratenum + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
                      polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
                      payments + ibrdcom + totcom, data = ibrdimp, 
                    index = c("ccode", "year"), model = "within",
                    effect = "twoways")
(seratefulllag1 <- coeftest(ratefulllag1, vcov=vcovHC(ratefulllag1, cluster="group"))) # interaction model

interaction_plot_continuous(model = ratefulllag1, conf = 0.9, effect = "nodefense", moderator = "ratenum", 
                            interaction = "nodefense:ratenum", title = "", xlabel = "Performance", 
                            ylabel = "Estimated Marginal Coefficient of No Alliance",
                            varcov = vcovHC(ratefulllag1, cluster="group")) #  exported and presented as Figure 4

ratefulllag2 <- plm(polrefben ~ defense + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
                      polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
                      payments + + ibrdcom + totcom, data = ibrdimp, 
                    index = c("ccode", "year"), model = "within",
                    effect = "twoways")
(seratefulllag2 <- coeftest(ratefulllag2, vcov=vcovHC(ratefulllag2, cluster="group"))) # ally only model

rpss <- plm(polrefben ~ nodefense + ratenum + numIOs + asymm + unsc +
              polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
              payments + ibrdcom + totcom, data = ibrdimp, 
            index = c("ccode", "year"), model = "within",
            effect = "twoways") # simple model without interaction for t test of interaction term
(serpss <- coeftest(rpss, vcov=vcovHC(rpss, cluster="group"))) # ally only model

waldtest(ratefulllag1, rpss, vcov = vcovHC) # p value is 0.011

ratefulllag4 <- plm(polrefben ~ defense + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
                      polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
                      payments + + ibrdcom + totcom, data = ibrdimp[!is.na(ibrdimp$ratenum),], 
                    index = c("ccode", "year"), model = "within",
                    effect = "twoways")
(seratefulllag4 <- coeftest(ratefulllag4, vcov=vcovHC(ratefulllag4, cluster="group"))) # just perf data

covs_full = c("Non-ally", "Performance", "Number IOs", "Vote-power asymmetry", "GDPPC",
           "GDP growth", "Population", "UNSC member", "Polity2", "UN voting (ideal pt dist from U.S.)", "U.S. aid", 
           "Trade with U.S.", "Reserves", "Receipts", "Payments", "World Bank commitments", "Bilateral and private commitments", 
           "Non-ally:Performance") # list of covariates for IBRD tests
covs_full2 = c("Ally", "Number IOs", "Vote-power asymmetry", "GDPPC",
              "GDP growth", "Population", "UNSC member", "Polity2", "UN voting (ideal pt dist from U.S.)", "U.S. aid", 
              "Trade with U.S.", "Reserves", "Receipts", "Payments", "World Bank commitments", "Bilateral and private commitments") # list of covariates for IBRD tests
ctryfe <- list(c("Country fixed effects", "Yes", "Yes", "Yes", "Yes")) # create rows denoting model specifications for regression tables
ctryre <- list(c("Country random effects", "Yes", "Yes", "Yes", "Yes"))
yearfe <- list(c("Year fixed effects", "Yes", "Yes", "Yes", "Yes"))
yearre <- list(c("Year random effects", "Yes", "Yes", "Yes", "Yes"))
wbfe <- list(c("World Bank President fixed effects", "Yes", "Yes", "Yes", "Yes"))
io <- list(c("Organization", "IBRD", "IBRD")) # create lists denoting relevant IOs

stargazer(ratebivlag1, ratepartlag1, ratefulllag1, 
          se = list(seratebivlag1[,2], seratepartlag1[,2], seratefulllag1[,2]),
          dep.var.labels = c("Number of reforms"),
          covariate.labels = covs_full,
          style = "ajps", omit = c("ccode", "year"), add.lines = c(ctryfe, yearfe, io),
          omit.stat = c("ll", "aic", "theta", "F", "rsq", "adj.rsq")) # Interaction results in table form; exported as Table 3

stargazer(ratebivlag2, ratepartlag2, ratefulllag2, 
          se = list(seratebivlag2[,2], seratepartlag2[,2], seratefulllag2[,2]),
          dep.var.labels = c("Number of reforms"),
          covariate.labels = covs_full2,
          style = "ajps", omit = c("ccode", "year"), add.lines = c(ctryfe, yearfe, io),
          omit.stat = c("ll", "aic", "theta", "F", "rsq", "adj.rsq")) # Ally results in table form; exported as Table 2

#### Robustness checks ####
##->## Country and year random effects
covs_rob <- c("Ally", "Non-ally", "Performance", "Number IOs", "Vote-power asymmetry", "GDPPC",
               "GDP growth", "Population", "UNSC member", "Polity2", "UN voting (ideal pt dist from U.S.)", "U.S. aid", 
               "Trade with U.S.", "Reserves", "Receipts", "Payments", "World Bank commitments", 
              "Bilateral and private commitments", "Non-ally:Performance") # list of covariates for robustness

stargazer(ratefulllag4, 
          se = list(seratefulllag4[,2]),
          dep.var.labels = c("Number of reforms"),
          covariate.labels = covs_full2,
          style = "ajps", omit = c("ccode", "year"), add.lines = c(ctryfe, yearfe, io),
          omit.stat = c("ll", "aic", "theta", "F", "rsq", "adj.rsq")) # Ally results for performance sample in table form; exported as Table A4, model 2

rob2 <- plm(polrefben ~ nodefense*ratenum + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
            polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
            payments + ibrdcom + totcom, data = ibrdimp, 
            index = c("ccode", "year"), model = "random", random.method = "amemiya",
            effect = "twoways")
(serob2 <- coeftest(rob2, vcov=vcovHC(rob2, cluster="group"))) # interaction model

wbrob2 <- plm(polrefben ~ defense + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
                polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
                payments + ibrdcom + totcom, data = ibrdimp, 
              index = c("ccode", "year"), model = "random", random.method = "amemiya",
              effect = "twoways")
(sewbrob2 <- coeftest(wbrob2, vcov=vcovHC(wbrob2, cluster="group"))) # ally model

stargazer(wbrob2, rob2, 
          se = list(sewbrob2[,2], serob2[,2]),
          dep.var.labels = c("Number of reforms"),
          covariate.labels = covs_rob,
          style = "ajps", omit = c("ccode", "year"), add.lines = c(ctryre, yearre, io),
          omit.stat = c("ll", "aic", "theta", "F", "rsq", "adj.rsq")) # exported as Table A4, models 4-5

##->## Alternate alliance measures

covs_alt <- c("Performance", "Number IOs", "Vote-power asymmetry", "GDPPC", "GDP growth", "Population", 
              "Polity2", "UN voting (ideal pt dist from U.S.)", "U.S. aid", "Trade with U.S.",
              "Reserves", "Receipts", "Payments", "World Bank commitments", 
              "Bilateral and private commitments", "Non-ally (entente):Performance",
              "Non-ally (non-aggression):Performance") # list alternate covariate labels
covs_alt2 <- c("Ally (entente)", "Ally (non-aggression)", "Non-ally (entente)", "Non-ally (non-aggression)",
              "Performance", "Number IOs", "Vote-power asymmetry", "GDPPC", "GDP growth", "Population", 
              "Polity2", "UN voting (ideal pt dist from U.S.)", "U.S. aid", "Trade with U.S.",
              "Reserves", "Receipts", "Payments", "World Bank commitments", 
              "Bilateral and private commitments", "Non-ally (entente):Performance",
              "Non-ally (non-aggression):Performance") # list alternate covariate labels
io3 <- list(c("Organization", "IBRD", "IBRD", "IBRD", "IBRD",
              "IBRD", "IBRD", "IBRD", "IBRD")) # additional IO list

rob3 <- plm(polrefben ~ noentente*ratenum + numIOs + asymm + gdppc + gdpgrow + pop + 
              polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + payments + ibrdcom + totcom, 
            data = ibrdimp, 
            index = c("ccode", "year"), model = "within",
            effect = "twoways")
(serob3 <- coeftest(rob3, vcov=vcovHC(rob3, cluster="group"))) # entente interaction

wbrob3 <- plm(polrefben ~ entente + numIOs + asymm + gdppc + gdpgrow + pop + 
              polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + payments + ibrdcom + totcom, 
              data = ibrdimp, 
            index = c("ccode", "year"), model = "within",
            effect = "twoways")
(sewbrob3 <- coeftest(wbrob3, vcov=vcovHC(wbrob3, cluster="group"))) # entente only

rob4 <- plm(polrefben ~ nononagg*ratenum + numIOs + asymm + gdppc + gdpgrow + pop + 
              polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + payments + ibrdcom + totcom, 
            data = ibrdimp, 
            index = c("ccode", "year"), model = "within",
            effect = "twoways")
(serob4 <- coeftest(rob4, vcov=vcovBK(rob4, cluster="group"))) # non-aggression interaction

wbrob4 <- plm(polrefben ~ nonaggression + numIOs + asymm + gdppc + gdpgrow + pop + 
              polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + payments + ibrdcom + totcom, 
              data = ibrdimp, 
            index = c("ccode", "year"), model = "within",
            effect = "twoways")
(sewbrob4 <- coeftest(wbrob4, vcov=vcovBK(wbrob4, cluster="group"))) # non-aggression only

stargazer(wbrob3, wbrob4, rob3, rob4, 
          se = list(sewbrob3[,2], sewbrob4[,2], serob3[,2], serob4[,2]),
          dep.var.labels = c("Number of reforms"),
          covariate.labels = covs_alt2,
          style = "ajps", omit = c("ccode", "year"), add.lines = c(ctryfe, yearfe, io3),
          omit.stat = c("ll", "aic", "theta", "F", "rsq", "adj.rsq")) # exported as Table A5

##->## No Imputation
covs_noimp <- c("Ally", "Non-ally", "Performance", "Number IOs", "Vote-power asymmetry", 
              "Population", "UNSC member", "Non-ally:Performance") # list of covariates for IBRD tests

rob9 <-  plm(polrefben ~ nodefense*ratenum + numIOs + asymm + pop + unsc, data = ibrdnoimp, 
             index = c("ccode", "year"), model = "within",
             effect = "twoways")       
(serob9 <- coeftest(rob9, vcov=vcovHC(rob9, cluster="group"))) # restricted covariates interaction

wbrob9 <-  plm(polrefben ~ defense + numIOs + asymm + pop + unsc, data = ibrdnoimp, 
             index = c("ccode", "year"), model = "within",
             effect = "twoways")       
(sewbrob9 <- coeftest(wbrob9, vcov=vcovHC(wbrob9, cluster="group"))) # restricted covariates ally

stargazer(wbrob9, rob9, 
          se = list(sewbrob9[,2], serob9[,2]),
          dep.var.labels = c("Number of reforms"),
          covariate.labels = covs_noimp,
          style = "ajps", omit = c("ccode", "year"), add.lines = c(ctryfe, yearfe, io),
          omit.stat = c("ll", "aic", "theta", "F", "rsq", "adj.rsq")) # exported as Table A6, models 5-6

##->## Only vote change reforms

rob10 <- lm(votechpos ~ nodefense*ratenum + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
              polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
              payments + ibrdcom + totcom + factor(year) + factor(ccode), data = ibrdimp)
(serob10 <- coeftest(rob10, vcov=vcovHC(rob10, type = "HC0", cluster=ibrdimp$ccode)))

rob10a <- lm(votechpos ~ defense + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
               polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
               payments + ibrdcom + totcom + factor(year) + factor(ccode), data = ibrdimp)
(serob10a <- coeftest(rob10a, vcov=vcovHC(rob10a, type = "HC0", cluster=ibrdimp$ccode)))

stargazer(rob10a, rob10, 
          se = list(serob10a[,2], serob10[,2]),
          dep.var.labels = c("Number of reforms"),
          covariate.labels = covs_rob,
          style = "ajps", omit = c("ccode", "year"), add.lines = c(ctryfe, yearfe, io),
          omit.stat = c("ll", "aic", "theta", "F", "rsq", "adj.rsq")) # noted in paper

##->## Logit model

cloglog = function(x) log(-log(1-x))

rob12 <- glm(polrefben ~ nodefense*ratenum + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
              polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
              payments + ibrdcom + totcom + factor(year) + factor(ccode), data = ibrdimp,
             family = binomial(link = cloglog), control=glm.control(maxit=200))
(serob12 <- coeftest(rob12, vcov=vcovHC(rob12, type = "HC0", cluster=ibrdimp$ccode)))

# the interaction model yields NA coefficients with fixed effects included. While the ally finding remains robust,
# we prefer OLS because it yields consistent results with coefficients for each variable. Similar
# problems manifest if the link is changed to "probit". Therefore, we use OLS, though we use GAM below.

rob12a <- glm(polrefben ~ defense + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
               polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
               payments + ibrdcom + totcom + factor(year) + factor(ccode), data = ibrdimp, 
              family = binomial(link = "logit"), control=glm.control(maxit=20))
(serob12a <- coeftest(rob12a, vcov=vcovHC(rob12a, type = "HC0", cluster=ibrdimp$ccode)))

stargazer(rob12a, 
          se = list(serob12a[,2]),
          dep.var.labels = c("Number of reforms"),
          covariate.labels = covs_full2,
          style = "ajps", omit = c("ccode", "year"), add.lines = c(ctryfe, yearfe, io),
          omit.stat = c("ll", "aic", "theta", "F", "rsq", "adj.rsq")) # Table A4, model 3

##->## Bank vs. borrower performance robustness check
iegfull <- read.csv("IEG_World_Bank_Project_Performance_Ratings.csv")
iegfull <- iegfull[c(5,19,28)] # subset data
iegfull$ccode <- countrycode(iegfull$Country.Name, origin = "country.name", destination = "cown") # create ccode
iegfull <- iegfull[-c(1)] # delete extra col
colnames(iegfull) <- c("year", "bankperf", "ccode") # change col names
iegfull$bankperf <- as.character(iegfull$bankperf) # change data format
iegfull <- iegfull[!(iegfull$bankperf == "NOT RATED" | 
                       iegfull$bankperf == "" | 
                       iegfull$bankperf == "NOT APPLICABLE" |
                       iegfull$bankperf == "NOT AVAILABLE"),] # subset to those projects for which Bank performance was rated
iegfull$bankperf <- ifelse(iegfull$bankperf == "HIGHLY SATISFACTORY", 6,
                           ifelse(iegfull$bankperf == "SATISFACTORY", 5, 
                                  ifelse(iegfull$bankperf == "MODERATELY SATISFACTORY", 4,
                                         ifelse(iegfull$bankperf == "MODERATELY UNSATISFACTORY", 3, 
                                                ifelse(iegfull$bankperf == "UNSATISFACTORY", 2, 1))))) # change coding to numeric
iegagg <- aggregate(iegfull$bankperf, by = list(iegfull$ccode, iegfull$year), FUN = mean) # take average perf by country-year
colnames(iegagg) <- c("ccode", "year", "bankperf")
ibrdimpi <- merge(unique(ibrdimp), unique(iegagg), by = c("ccode", "year")) # merge data

sum(ibrdimpi$polrefben) # see how many reforms are left in the subsetted data

# Here, we try to use the IEG rating that measures purely Bank performance. Only 3 cases of beneficial 
# reforms remain from the broader sample, so we are unable to use it. We note this in text.

##->## WB Pres FEs
wbpres <- read.csv("wb_pres.csv")
ibrdimpj <- merge(unique(ibrdimp), unique(wbpres), by = c("year")) # merge data
colnames(ibrdimpj)[74] <- "wbpresident"

rob14 <- lm(polrefben ~ nodefense*ratenum + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
               polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
               payments + ibrdcom + totcom + factor(wbpresident) + factor(ccode), data = ibrdimpj)
(serob14 <- coeftest(rob14, vcov=vcovHC(rob14, type = "HC0", cluster=ibrdimpj$ccode))) # interaction model

rob14a <- lm(polrefben ~ defense + numIOs + asymm + gdppc + gdpgrow + pop + unsc + polity2 + 
                absidealimportantdiff + USaid + tottrade + reserves + receipts + 
                payments + ibrdcom + totcom + factor(wbpresident) + factor(ccode), data = ibrdimpj)
(serob14a <- coeftest(rob14a, vcov=vcovHC(rob14a, type = "HC0", cluster=ibrdimpj$ccode))) # ally only model

stargazer(rob14a, rob14, 
          se = list(serob14a[,2], serob14[,2]),
          dep.var.labels = c("Number of reforms"),
          covariate.labels = covs_rob,
          style = "ajps", omit = c("ccode", "wbpresident"), add.lines = c(ctryfe, wbfe, io),
          omit.stat = c("ll", "aic", "theta", "F", "rsq", "adj.rsq")) # Table A6, models 7-8

##->## Malik and Stone performance

malikstone <- read.csv("Malik_Stone_Full.csv") # load data from Malik and Stone for alternate performance measure
malikstone <- na.omit(malikstone) # eliminate NAs for aggregation
malikstone$ccode <- countrycode(malikstone$ctry, origin = "country.name", destination = "cown")
malikstone <- aggregate(list(malikstone$avgPDO), by = list(malikstone$year, malikstone$ccode), FUN = mean)
colnames(malikstone) <- c("year", "ccode", "msgrade")
malikstone <- merge(unique(malikstone), ibrdimp, by = c("ccode", "year"))

# Here, we try to use the alternate performance measure from Malik and Stone (2016). 0 cases of beneficial 
# reforms remain from the broader sample, so we are unable to use it. We note this in text.

##->## Lag DV 

ibrdimpk <- 
  ibrdimp %>%
  group_by(ccode) %>%
  mutate(lag.polrefben = dplyr::lag(polrefben, n = 1, default = NA))

rob15 <- lm(polrefben ~ nodefense*ratenum + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
              polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
              payments + ibrdcom + totcom + lag.polrefben + factor(year) + factor(ccode), data = ibrdimpk)
(serob15 <- coeftest(rob15, vcov=vcovHC(rob15, type = "HC0", cluster=ibrdimpk$ccode)))

rob15a <- lm(polrefben ~ defense + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
               polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
               payments + ibrdcom + totcom + lag.polrefben + factor(year) + factor(ccode), data = ibrdimpk)
(serob15a <- coeftest(rob15a, vcov=vcovHC(rob15a, type = "HC0", cluster=ibrdimpk$ccode)))

covs_alt3 <- c("Ally", "Non-ally", "Performance", "Number IOs", "Vote-power asymmetry", "GDPPC",
              "GDP growth", "Population", "UNSC member", "Polity2", "UN voting (ideal pt dist from U.S.)", "U.S. aid", 
              "Trade with U.S.", "Reserves", "Receipts", "Payments", "World Bank commitments", 
              "Bilateral and private commitments", "Lag DV", "Non-ally:Performance") # list of covariates

stargazer(rob15a, rob15, 
          se = list(serob15a[,2], serob15[,2]),
          dep.var.labels = c("Number of reforms"),
          covariate.labels = covs_alt3,
          style = "ajps", omit = c("ccode", "year"), add.lines = c(ctryfe, yearfe, io),
          omit.stat = c("ll", "aic", "theta", "F", "rsq", "adj.rsq")) # Table A6, models 1-2

##->## Drop Reforms US can veto (vote share)

ibrdimpv <- ibrdimp # copy data frame
ibrdimpv$polrefben2 <- ibrdimpv$polrefben - ibrdimpv$votechpos # subtract vote change from reforms
ibrdimpv$polrefben2[ibrdimpv$polrefben2 < 1] <- 0 # change negatives to 0

rob16 <- lm(polrefben2 ~ nodefense*ratenum + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
              polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
              payments + ibrdcom + totcom + factor(year) + factor(ccode), data = ibrdimpv)
(serob16 <- coeftest(rob16, vcov=vcovHC(rob16, type = "HC0", cluster=ibrdimpv$ccode)))

rob16a <- lm(polrefben2 ~ defense + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
               polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
               payments + ibrdcom + totcom + factor(year) + factor(ccode), data = ibrdimpv)
(serob16a <- coeftest(rob16a, vcov=vcovHC(rob16a, type = "HC0", cluster=ibrdimpv$ccode)))

stargazer(rob16a, rob16, 
          se = list(serob16a[,2], serob16[,2]),
          dep.var.labels = c("Number of reforms"),
          covariate.labels = covs_rob,
          style = "ajps", omit = c("ccode", "year"), add.lines = c(ctryfe, yearfe, io),
          omit.stat = c("ll", "aic", "theta", "F", "rsq", "adj.rsq")) # Table A6, models 9-10

##->## Rolling average on performance 

ibrdimp2 <- ibrdimp
for (i in 1:7845) {
  yrtemp <- ibrdimp2$year[i]
  temp <- mean(ibrdimp2$ratenum[ibrdimp2$year == yrtemp], na.rm=T)
  ibrdimp2$ratenum[i] <- ifelse(is.na(ibrdimp2$ratenum[i]), temp, ibrdimp2$ratenum[i])
} # replace missing performance values with year average
ibrdimp2$rollperf <- rollmean(ibrdimp2$ratenum, k = 5, fill=NA) # compute five-year rolling average

rob17 <- lm(polrefben ~ nodefense*rollperf + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
              polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
              payments + ibrdcom + totcom + factor(year) + factor(ccode), data = ibrdimp2)
(serob17 <- coeftest(rob17, vcov=vcovHC(rob17, type = "HC0", cluster=ibrdimp2$ccode)))

stargazer(rob17, 
          se = list(serob17[,2]),
          dep.var.labels = c("Number of reforms"),
          covariate.labels = covs_full,
          style = "ajps", omit = c("ccode", "year"), add.lines = c(ctryfe, yearfe, io),
          omit.stat = c("ll", "aic", "theta", "F", "rsq", "adj.rsq")) # Table A4, model 1

##->## Collapse IEG ratings to 0,1 scale

ibrdimp3 <- ibrdimp
ibrdimp3$ratenumbin <- ifelse(ibrdimp3$ratenum > 3, 1, 0)

rob25 <- plm(polrefben ~ nodefense*ratenumbin + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
               polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
               payments + ibrdcom + totcom, data = ibrdimp3, 
             index = c("ccode", "year"), model = "within",
             effect = "twoways")
(serob25 <- coeftest(rob25, vcov=vcovHC(rob25, cluster="group"))) # interaction model; Table A6, model 11

##->## Eliminate watershed reforms
# 1959: IDA creation
# 1985: MIGA creation
# 1994: environment and sustainable development shifts
# 2010-11: voice and participation reforms

ibrdimp3 <- ibrdimp
ibrdimp3 <- ibrdimp3[!(ibrdimp3$year %in% c(1960, 1986, 1995, 2011, 2012)),] # excluding watershed

rob18 <- plm(polrefben ~ nodefense*ratenum + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
                      polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
                      payments + ibrdcom + totcom, data = ibrdimp3, 
                    index = c("ccode", "year"), model = "within",
                    effect = "twoways")
(serob18 <- coeftest(rob18, vcov=vcovHC(rob18, cluster="group"))) # interaction model

rob18a <- plm(polrefben ~ defense + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
                      polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
                      payments + + ibrdcom + totcom, data = ibrdimp3, 
                    index = c("ccode", "year"), model = "within",
                    effect = "twoways")
(serob18a <- coeftest(rob18a, vcov=vcovHC(rob18a, cluster="group"))) # ally only model

stargazer(rob18a, rob18, 
          se = list(serob18a[,2], serob18[,2]),
          dep.var.labels = c("Number of reforms"),
          covariate.labels = covs_rob,
          style = "ajps", omit = c("ccode", "year"), add.lines = c(ctryfe, yearfe, io),
          omit.stat = c("ll", "aic", "theta", "F", "rsq", "adj.rsq")) # Table A6, models 3-4

## Other measures of alignment ##

## IO membership overlap
joint <- read.csv("dyadic_formatv3.csv") # load data from COW (v3) on dyadic state membership in IOs
joint <- joint[joint$ccode1==2,] # restrict to US overlap
joint[joint==-1] <- 0 # change all 0 to 0
joint[joint==-9] <- 0 # change all 0 to 0

state <- read.csv("state_year_formatv3.csv")
state <- state[state$ccode == 2,]
state[state==-1] <- 0 # change all 0 to 0
state[state==-9] <- 0 # change all 0 to 0
state$numjoined <- rowSums(state[ , c(4:537)])
state <- state[c(2,538)]

joint$ioverlap <- rowSums(joint[ , c(7:540)])
joint <- joint[c(3,5,541)]
colnames(joint)[1] <- "ccode"
joint <- left_join(joint, state, by = c("year"))
joint$year <- joint$year+1
joint$shareoverlap <- joint$ioverlap / joint$numjoined
joint <- left_join(unique(ibrdimp), unique(joint), by = c("ccode", "year"))
joint$noverlap <- 1 - joint$shareoverlap

rob2p <- plm(polrefben ~ noverlap*ratenum + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
               polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
               payments + ibrdcom + totcom, data = joint, 
             index = c("ccode", "year"), model = "within",
             effect = "twoways")
(serob2p <- coeftest(rob2p, vcov=vcovHC(rob2p, cluster="group"))) # interaction model

joint$shareoverlap[is.na(joint$shareoverlap)] <- 0
wbrob2p <- plm(polrefben ~ shareoverlap + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
                 polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
                 payments + ibrdcom + totcom, data = joint, 
               index = c("ccode", "year"), model = "within",
               effect = "twoways")
(sewbrob2p <- coeftest(wbrob2p, vcov=vcovHC(wbrob2p, cluster="group"))) # ally model
## Also in Table A5

covs_rob <- c("IO overlap", "IO non-overlap", "Performance", "Number IOs", "Vote-power asymmetry", "GDPPC",
              "GDP growth", "Population", "UNSC member", "Polity2", "UN voting (ideal pt dist from U.S.)", "U.S. aid", 
              "Trade with U.S.", "Reserves", "Receipts", "Payments", "World Bank commitments", 
              "Bilateral and private commitments", "Non-ally:Performance") # list of covariates for robustness
 
## Dyadic trade with US ##
trade <- read.csv("Dyadic_COW_4.0.csv") # load data v4.0 from cow on dyadic trade
trade <- trade[trade$ccode1==2,] # subet to US only
trade <- trade[-c(1)] # delete first column
colnames(trade)[1] <- "ccode"
trade <- trade[c(1,2,9)]
trade$year <- trade$year + 1
 
trade <- left_join(unique(ibrdimp), unique(trade), by = c("ccode", "year"))
trade$smoothtotrade <- log1p(trade$smoothtotrade)
trade$smoothtotrade[is.na(trade$smoothtotrade)] <- 0
 
rob2t <- plm(polrefben ~ smoothtotrade*ratenum + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
                polity2 + absidealimportantdiff + USaid + reserves + receipts + 
                payments + ibrdcom + totcom, data = trade, 
              index = c("ccode", "year"), model = "within",
              effect = "twoways")
(serob2t <- coeftest(rob2t, vcov=vcovHC(rob2t, cluster="group"))) # interaction model
 
wbrob2t <- plm(polrefben ~ smoothtotrade + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
                  polity2 + absidealimportantdiff + USaid + reserves + receipts + 
                  payments + ibrdcom + totcom, data = trade, 
                index = c("ccode", "year"), model = "within",
                effect = "twoways")
(sewbrob2t <- coeftest(wbrob2t, vcov=vcovHC(wbrob2t, cluster="group"))) # ally model

covs_alt2 <- c("Ally (entente)", "Ally (non-aggression)", "IO member overlap", "Trade with U.S.",
               "Non-ally (entente)", "Non-ally (non-aggression)", "IO member non-overlap",
               "Performance", "Number IOs", "Vote-power asymmetry", "GDPPC", "GDP growth", "Population", "UNSC member",
               "Polity2", "UN voting (ideal pt dist from U.S.)", "U.S. aid", "Trade with U.S.",
               "Reserves", "Receipts", "Payments", "World Bank commitments", 
               "Bilateral and private commitments", "Non-ally (entente) X Performance",
               "Non-ally (non-aggression) X Performance", "IO member non-overlap X Performance",
               "Trade with U.S. X Performance") # list alternate covariate labels

stargazer(wbrob3, wbrob4, wbrob2p, wbrob2t, 
          rob3, rob4, rob2p, rob2t,  
          se = list(sewbrob3[,2], sewbrob4[,2], sewbrob2p[,2], sewbrob2t[,2], 
                    serob3[,2], serob4[,2], serob2p[,2], serob2t[,2]),
          dep.var.labels = c("Number of reforms"),
          covariate.labels = covs_alt2,
          style = "ajps", omit = c("ccode", "year"), add.lines = c(ctryfe, yearfe, io3),
          omit.stat = c("ll", "aic", "theta", "F", "rsq", "adj.rsq")) # exported as Table A5; combine with other measures of alignment for appendix table
 
## Diffusion controls
# Time trends
ibrdimp$year2 <- ibrdimp$year^2
ibrdimp$year3 <- ibrdimp$year^3
rob2d <- lm(polrefben ~ nodefense*ratenum + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
               polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
               payments + ibrdcom + totcom + year + year2 + year3 + factor(ccode), data = ibrdimp)
(serob2d <- coeftest(rob2d, vcov=vcovHC(rob2d, type="HC0", cluster=ibrdimp$ccode))) # interaction model

wbrob2d <- lm(polrefben ~ defense + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
                 polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
                 payments + ibrdcom + totcom + year + year2 + year3 + factor(ccode), data = ibrdimp)
(sewbrob2d <- coeftest(wbrob2d, vcov=vcovHC(wbrob2d, type="HC0", cluster=ibrdimp$ccode))) # ally model

# Avg reforms by category
cats <- read.csv("reformp_cat_ibrd.csv")
cats <- aggregate(list(cats$vote, cats$miss, cats$rule, cats$staff, cats$lend), by = list(cats$year), FUN = mean)
colnames(cats) <- c("year", "vote", "miss", "rule", "staff", "lend")
ibrdimp <- left_join(ibrdimp, cats, by = "year")
rob2c <- plm(polrefben ~ nodefense*ratenum + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
               polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
               payments + ibrdcom + totcom + vote + miss + rule + staff + lend, data = ibrdimp, 
             index = c("ccode"), model = "within")
(serob2c <- coeftest(rob2c, vcov=vcovHC(rob2c, cluster="group"))) # interaction model

wbrob2c <- plm(polrefben ~ defense + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
                 polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
                 payments + ibrdcom + totcom + vote + miss + rule + staff + lend, data = ibrdimp, 
               index = c("ccode"), model = "within")
(sewbrob2c <- coeftest(wbrob2c, vcov=vcovHC(wbrob2c, cluster="group"))) # ally model

covs_alt3 <- c("Ally", "Performance", "Non-ally",
               "Performance", "Number IOs", "Vote-power asymmetry", "GDPPC", "GDP growth", "Population", "UNSC member",
               "Polity2", "UN voting (ideal pt dist from U.S.)", "U.S. aid", "Trade with U.S.",
               "Reserves", "Receipts", "Payments", "World Bank commitments", 
               "Bilateral and private commitments", "Non-ally X Performance",
               "Non-ally (non-aggression) X Performance", "IO member non-overlap X Performance",
               "Trade with U.S. X Performance") # list alternate covariate labels

iocon <- list(c("Average # reforms by category and year","Yes", "No", "Yes", "No"))
tt <- list(c("Quadratic time trend","No", "Yes", "No", "Yes"))
stargazer(wbrob2c, wbrob2d, rob2c, rob2d,
          se = list(sewbrob2c[,2], sewbrob2d[,2], serob2c[,2], serob2d[,2]),
          dep.var.labels = c("Number of reforms"),
          covariate.labels = covs_rob,
          style = "ajps", omit = c("ccode", "year", "year", "year2", "year3",
                                   "vote", "miss", "rule", "staff", "lend"), add.lines = c(ctryfe, iocon, tt),
          omit.stat = c("ll", "aic", "theta", "F", "rsq", "adj.rsq")) # exported as Table A7

##->## Adjust units
unitsalt <- ibrdimp
unitsalt$numIOs <- scale(unitsalt$numIOs)
unitsalt$asymm <- scale(unitsalt$asymm)
unitsalt$gdppc <- scale(unitsalt$gdppc)
unitsalt$pop <- scale(unitsalt$pop)
unitsalt$gdpgrow <- scale(unitsalt$gdpgrow)
unitsalt$polity2 <- scale(unitsalt$polity2)
unitsalt$absidealimportantdiff <- scale(unitsalt$absidealimportantdiff)
unitsalt$USaid <- scale(unitsalt$USaid)
unitsalt$totcom <- scale(unitsalt$totcom)
unitsalt$ibrdcom <- scale(unitsalt$ibrdcom)
unitsalt$reserves <- scale(unitsalt$reserves)
unitsalt$receipts <- scale(unitsalt$receipts)
unitsalt$payments <- scale(unitsalt$payments)
unitsalt$ratenum <- scale(unitsalt$ratenum)

robun <-  plm(polrefben ~ defense + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
                polity2 + absidealimportantdiff + USaid + tottrade + reserves + 
                receipts + payments + ibrdcom + totcom, 
              data = unitsalt, 
              index = c("ccode", "year"), model = "within",
              effect = "twoways")       
(serobun <- coeftest(robun, vcov=vcovHC(robun, cluster="group")))

robuna <-  plm(polrefben ~ nodefense*ratenum + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
               polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + payments + ibrdcom + totcom, 
             data = unitsalt, 
             index = c("ccode", "year"), model = "within",
             effect = "twoways")       
(serobuna <- coeftest(robuna, vcov=vcovHC(robuna, cluster="group")))

stargazer(robun, robuna, 
          se = list(serobun[,2], serobuna[,2]),
          dep.var.labels = c("Number of reforms"),
          covariate.labels = covs_rob,
          style = "ajps", omit = c("ccode", "year"), add.lines = c(ctryfe, yearfe, io),
          omit.stat = c("ll", "aic", "theta", "F", "rsq", "adj.rsq")) # Table A8

##->## GDP split sample
quantile(ibrdimp$gdp)
ibrdimpl <- ibrdimp[ibrdimp$gdp < 11.66,] # drop lowest quartile
robgdpl <- plm(polrefben ~ nodefense*ratenum + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
                      polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
                      payments + ibrdcom + totcom, data = ibrdimpl, 
                    index = c("ccode", "year"), model = "within",
                    effect = "twoways")
(segdpl <- coeftest(robgdpl, vcov=vcovHC(robgdpl, cluster="group"))) # interaction model

robgdpl2 <- plm(polrefben ~ defense + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
                      polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
                      payments + + ibrdcom + totcom, data = ibrdimpl, 
                    index = c("ccode", "year"), model = "within",
                    effect = "twoways")
(segdpl2 <- coeftest(robgdpl2, vcov=vcovHC(robgdpl2, cluster="group"))) # ally only model

stargazer(robgdpl2, robgdpl,
          se = list(segdpl2[,2], segdpl[,2]),
          dep.var.labels = c("Number of reforms"),
          covariate.labels = covs_rob,
          style = "ajps", omit = c("ccode", "year"), add.lines = c(ctryfe, yearfe, io),
          omit.stat = c("ll", "aic", "theta", "F", "rsq", "adj.rsq")) # Table A9

##->## Non-parametric ##
int.gam.all <- gam(polrefben ~ s(ratenum, bs = "cr", by = nodefense) + numIOs + asymm + gdppc + gdpgrow + pop + unsc +
                     polity2 + absidealimportantdiff + USaid + tottrade + reserves + receipts + 
                     payments + ibrdcom + totcom,
               family = binomial, data = ibrdimp)
summary(int.gam.all)

int.gam <- gam(polrefben ~ s(ratenum, bs = "cr", by = nodefense) + asymm + gdppc +
                 polity2 + absidealimportantdiff, 
                 family = binomial, data = ibrdimp)
summary(int.gam)

# We limit the covariates to those that are substantively most important 
# and show significance most often in our main models, though results for 
# the interaction term are very similar with the full cohort. We limit them 
# in order to construct plots below; we run into memory issues plotting otherwise.

xrng <- list(ratenum = range(ibrdimp$ratenum, na.rm = TRUE))
ratenum.seq <- seq(xrng$ratenum[1], xrng$ratenum[2], length.out = 25)
nodefense.seq <- c(1,0)
xrng <- list(gdppc = range(ibrdimp$gdppc, na.rm = TRUE))
gdppc.seq <- seq(xrng$gdppc[1], xrng$gdppc[2], length.out = 25)
xrng <- list(asymm = range(ibrdimp$asymm, na.rm = TRUE))
asymm.seq <- seq(xrng$asymm[1], xrng$asymm[2], length.out = 25)
xrng <- list(polity2 = range(ibrdimp$polity2, na.rm = TRUE))
polity2.seq <- seq(xrng$polity2[1], xrng$polity2[2], length.out = 25)
xrng <- list(absidealimportantdiff = range(ibrdimp$absidealimportantdiff, na.rm = TRUE))
absidealimportantdiff.seq <- seq(xrng$absidealimportantdiff[1], xrng$absidealimportantdiff[2], length.out = 25)
xrng <- list(USaid = range(ibrdimp$USaid, na.rm = TRUE))

newdata <- expand.grid(nodefense = nodefense.seq, 
                       asymm = asymm.seq,
                       ratenum = ratenum.seq,
                       gdppc = gdppc.seq, 
                       polity2 = polity2.seq,
                       absidealimportantdiff = absidealimportantdiff.seq)
yhat.int <- predict(int.gam, newdata, type = "response")

plot(ratenum.seq, yhat.int[1:25], ylim = c(0, 0.5), t = "l", col = "red",
     xlab = "World Bank Performance", ylab = "Probability of Benefiting from Reform",
     main = "")
lines(ratenum.seq, yhat.int[1:25 + 25], lty = 2, col = "blue", type = "s")
legend("bottomleft", legend = c("Non-ally", "Ally"), lty = 1:2,
       col = c("red", "blue"), horiz = TRUE, bty = "n")

pdata <- data.frame(ratenum = ratenum.seq, yhat = yhat.int[1:25],
                    type = c(rep("Non-ally", 25)))

# Raw data
ggplot(pdata, aes(x = ratenum, y = yhat)) +
  geom_line(stat = "smooth") +
  ylab("Prob Non-ally Benefits from Reform") +
  xlab("IEG Performance") +
  theme(axis.title.x = element_text(size = 22), axis.title.y = element_text(size = 22),
        axis.text.x = element_text(size = 20), axis.text.y = element_text(size = 20),
        legend.text = element_text(size = 20), legend.title = element_text(size = 22))
## Table Figure A4  

#### Generalizability ####
dp <- read.csv("SIGO.csv") # load replication file from Davis and Pratt (2020) in ROIO
dp <- unique(dp[c(5,7:11)])
dp <- dp[sp$Security_IGO==0,]
dp <- dp[-c(3)] # exclude security IOs
dp_sub <- sample(dp$IGO, 100)
dp_sub <- dp[dp$IGO %in% dp_sub,]
cowio <- read_dta("igo_year_format_3.dta") # Load v3 of COW IGO data
cowio <- cowio[c(1:2)]
colnames(cowio) <- c("IGO", "IGO_Name")
dp_sub <- left_join(dp_sub, unique(cowio), by = c("IGO"))
# write.csv(dp_sub, "davis_pratt_io_list.csv") # write IO list to file; hand-coded our four
# generalizability categories for it
dp_full <- read.csv("davis_pratt_io_list.csv") # load data with hand coding
dp_full$general <- ifelse(dp_full$weighted_voting == 1 & dp_full$mps == 1 
                          & dp_full$performance == 1 & dp_full$exit == 1, 1, 0)
mean(dp_full$general) # 17%
mean(dp_full$general[dp_full$Econ_IGO == 1]) # 21% for econ IOs
mean(dp_full$general[dp_full$Social_IGO == 1]) # 5% for econ IOs
mean(dp_full$general[dp_full$Environ_IGO == 1]) # 18% for environmental IOs
mean(dp_full$general[dp_full$regional_org == 1]) # 22% for regional orgs

mean(dp_full$weighted_voting) # 19%
mean(dp_full$mps) # 28%
mean(dp_full$performance) # 39%
mean(dp_full$exit) # 31%

mean(dp_full$weighted_voting[dp_full$Econ_IGO == 1]) # 22%
mean(dp_full$mps[dp_full$Econ_IGO == 1]) # 33%
mean(dp_full$performance[dp_full$Econ_IGO == 1]) # 46%
mean(dp_full$exit[dp_full$Econ_IGO == 1]) # 38%

mean(dp_full$weighted_voting[dp_full$Environ_IGO == 1]) # 18%
mean(dp_full$mps[dp_full$Environ_IGO == 1]) # 27%
mean(dp_full$performance[dp_full$Environ_IGO == 1]) # 36%
mean(dp_full$exit[dp_full$Environ_IGO == 1]) # 18%

mean(dp_full$weighted_voting[dp_full$Social_IGO == 1]) # 5%
mean(dp_full$mps[dp_full$Social_IGO == 1]) # 18%
mean(dp_full$performance[dp_full$Social_IGO == 1]) # 18%
mean(dp_full$exit[dp_full$Social_IGO == 1]) # 9%

mean(dp_full$weighted_voting[dp_full$regional_org == 1]) # 24%
mean(dp_full$mps[dp_full$regional_org == 1]) # 29%
mean(dp_full$performance[dp_full$regional_org == 1]) # 44%
mean(dp_full$exit[dp_full$regional_org == 1]) # 41%

### These figures appear in Table 4

